Critical Behavior in the Gravitational Collapse of a Scalar Field with Angular 

Momentum in Spherical Symmetry 
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We study the critical collapse of a massless scalar field with angular momentum in spherical 
symmetry. In order to mimic the effects of angular momentum we perform a sum of the stress- 
energy tensors for all the scalar fields with the same eigenvalue I of the angular momentum operator 
and calculate the equations of motion for the radial part of these scalar fields. We have found 
that the critical solutions for different values of I are discretely self-similar (as in the original / = 
case). The value of the discrete, self-similar period, Aj, decreases as I increases in such a way that 
the critical solution appears to become periodic in the limit. The mass scaling exponent, 71, also 
decreases with I. 



I. INTRODUCTION 



Most studies of black hole critical phenomena (see [l[ , 
for reviews) to date (or related phenomena in other 
sets of nonlinear evolution equations) have been per- 
formed assuming spherical symmetry as a simplifyingas- 
sumption (exceptions are [3], [1] and more recently [3]) ■ 
This simplification has been adopted in most cases be- 
cause accurate calculation of Type II critical solutions — 
which exhibit structure at all scales due to their self- 
similar nature — requires great computational resources. 
Since spherically symmetric spacetimes do not allow for 
angular momentum, very little is currently known about 
the role of angular momentum in critical collapse. For 
a few cases, most notably the Type II solutions found 
in spherically symmetric collapse of a massless scalar 
field [S| , or certain types of perfect fluid |3] , [l^ , pertur- 
bative calculations about the spherical critical solutions 
suggest that non-spherical modes, including those con- 
tributing to net angular momentum, arc damped as one 
approaches criticality-'^ . In particular in j^, Q and 
using second order perturbation theory it was predicted 
that the angular momentum of the black holes produced 
should have the following dependence as a function of the 
critical parameter p: 



(1) 
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^ There is some numerical evidence for growing non-spherical 
modes in near-critical collapse, both for perfect fluids [lOll and 
massless scalar fields Q. However, these modes appear to grow 
so slowly that, in both cases, it is expected that the spherical 
unstable mode continues to dominate near criticality. 



where Lq is family-dependent and /x is a universal scaling 
exponent satisfying ^ > 27 (7 being the scaling exponent 
for the black hole mass). Specifically, it was suggested 
that jjL « 0.76 for the scalar field case, whereas the com- 
putations indicated that /i would depend on the equation 
of state for perfect fluid collapse. These calculations thus 
suggest that, at least for small deviations from spherical 
symmetry, the resulting solutions on the verge of black 
hole formation should remain spherically symmetric in 
non-symmetric collapse. We also note that an axisym- 
metric numerical relativity code has been developed 
to study non-perturbatively some effects of angular mo- 
mentum in the critical collapse of a scalar field. Inter- 
estingly, the results found for A and 7 in the case of a 
complex scalar field with principal azimuthal "quantum 
number" , m — 1 are very close to the results we find in 
our model for / = 1, as described in Sec. IIIII 



Here a different approach is taken. Maintaining spher- 
ical symmetry, the equations of motion for a massless 
scalar field are modified by effective terms which mock 
up some of the effects of angular momentum. As de- 
scribed below, the procedure amounts to performing an 
angular average over the matter field variables — similar 
to that done in [l^], and 14 1 — and results in an entire 
family of models, parameterized by a principal angular 
"quantum number" , / (we will generally restrict I to take 
on non-negative integer values, although real- valued Vs 
are also formally possible). We note that since the mod- 
els remain spherically symmetric, we cannot use them 
to address the validity of the perturbative calculations 
mentioned above (e.g. equation ([1])). Nonetheless, we 
find interesting results that may shed some light on the 
effects of angular momentum near the black hole thresh- 
old. 

Some of the main results that have been found are as 
follows. First, each value of the angular momentum pa- 
rameter I apparently defines a distinct critical solution. 
For / < 10, these solutions are found to be discretely self 
similar, with values of the echoing exponent. A;, that 
rapidly decrease (approximately exponentially) as I in- 
creases. As a result, for large values of and for the 
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time scales for which we are able to dynamically evolve 
near criticality, the threshold solutions become approxi- 
mately periodic. In addition, and as expected for Type 
II solutions, we find that for Z < 7 the masses of the 
black holes formed follow power laws. As with the echo- 
ing exponents, for increasing values of I it is found that 
the mass-scaling exponent, 7;, rapidly decreases, again 
approximately exponentially in I. 

The remainder of this paper is structured as follows. 
In the following section we describe the recipe used to 
calculate the effective equations of motion, along with 
the regularity and boundary conditions imposed in the 
solution of these equations. In Sec. IIIII we briefly de- 
scribe the numerical code, the way the solutions have 
been analyzed, and then provide a summary of the re- 
sults obtained for varying values of I. Throughout this 
paper we use units such that the universal gravitational 
constant, G, and the speed of light in vacuum, c, are both 
unity. 



where gab is the metric of the spacetime and Va is the 
metric-compatible covariant derivative. Again by con- 
struction, and as is proven in Appendix |Al the sum of 
these stress tensors 



ahi 
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is spherically symmetric, and thus depends only on 
ijj{t,r), I, and the metric Qab- We can now compute the 
effective equation of motion for the field, ip(t,r), using 
the fact that the divergence of the total stress energy 
tensor is zero, as is also proven in Appendix [Xl 



(6) 



The equations for the geometric variables are determined 
from the 3-1-1 decomposition of the Einstein field equa- 
tions. For the current study we adopt Schwarzschild-like 
(polar- areal) coordinates, in which the metric takes the 
form: 



II. EQUATIONS OF MOTION 



ds^ ^ -a^{t,r)dt^ + a^{t,r)dr'^ 



r^de^ 



2 • 2 

r sm 



0d(j)^ . (7) 



A. Equations 

In order to derive equations of motion, scalar fields of 
the following form are considered: 



m — —I, ~l + I — 1,1, 



(2) 



where Qim{0, <j)) are normalized rea/ eigenfunctions of the 
angular part of the fiatspace Laplacian with eigenvalue 
l{l + 1), and the index m labels the 2Z -I- 1 distinct or- 
thonormal eigcnfunctions for a given value of More 
explicitly: 
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TO = 0, 
TO > 0, 

TO < 0, 



(3) 

where Yim {0, </>) are the regular spherical harmonics. By 
construction, the scalar fields are not, in general, 
spherically symmetric and we therefore do not study 
their collapse directly. Instead, our strategy is to find ef- 
fective equations for the single {t, r)-dependent quantity 
'(/'^'•'(i, r), which we hereafter denote simply by ijj. To do 
so, for a specific value of we consider the stress-energy 
tensors for the 21 + 1 fields ^-I": 



ah 



^<?a6(v^*rvc^r), (4) 



^ Note that, in general, Qim{d, fp) will not be eigenfunctions of the 
azimuthal rotation operator {d/d4>) since they are real. 



Here a{t, r) is the lapse function and a{t, r) is the only 
non-trivial components of the 3-metric (both a and a 
are positive functions). Using this metric, the non zero 
components of the stress-energy tensor for a general value 
of I are 
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(n2 + $2) + 
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{21 + 1) 



'^')-l{l + l)^ 



and the stress-energy trace is 
{21 + 1) 



,(8) 
(9) 
(10) 

(11) 
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n2-$2j ^2l{l + l) 



21 



,(12) 



In the above expressions, we have made use of the auxil- 
iary variables, and 11, defined as follows: 
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(13) 
(14) 



The dynamical equations of motion for these fields, which 
follow from the definition of $ as well as the wave equa- 
tion for -0 (which in turn can be derived from the van- 
ishing of the divergence of the total stress tensor ^) are 
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then: 
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dr^a ) ' 
r"^ dr \ a ' 



l{l + l)aa 



(15) 
(16) 



Note that the dependence of these equations on I is 
only through the last term in equation (|16p which is 
proportional to Z(/ + l)/?'^- This term can be thought 
of as the field-theoretic extension of an analogous term 
due to the angular momentum potential, P /r^, in the 
1-dimensional reduced problem of a particle moving in a 
central potential. 

As mentioned above, equations for the geometric vari- 
ables result from the 3-1-1 decomposition of the field equa- 
tions, as well as from our choice of coordinates. Specifi- 
cally, we have the following 
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r n^ + $^ + /(/ + i)— V-^ 



r n2 + a>2-/(/ + i)^V' 



2r ' 
{21 + l)ran$. 



(17) 

(18) 
(19) 



Equation (|17p is the Hamiltonian constraint, which is 
used to determine the 3-metric component, a. Similarly, 
the slicing condition fixes the lapse function a at 
each instant of time, and is often known as the polar slic- 
ing condition. It can be derived from the demand that 
Tr [Kab) = K\ + K\ + K-^^ = K\ + 2K\ = 0, for 
all times. The Hamiltonian constraint and slicing con- 
dition, with appropriate regularity and boundary con- 
ditions, completely fix the geometric variables in this 
coordinate system. Equation (|19p is an extra equation 
derived from the definition of r and the momentum 
constraint. In our numerical solutions, it is used as a 
gauge of the accuracy of our calculations, as well as to 
provide a replacement for the Hamiltonian constraint in 
certain strong field instances where the numerical con- 
straint solver fails. In addition, we compute the mass 
aspect function, M(t,r), 



(20) 



which serves as a valuable diagnostic quantity in our cal- 
culations. The value of this function as r ^ oo agrees 
with the ADM mass, and more generally, in a vacuum 
region of spacetime, measures the amount of (gravitat- 
ing) mass contained within the 2-sphere of radius r at 
time t. Moreover, 2M{t,r)/r is useful since its value 
approaches 1 when a trapped surface is developing and 



hence (modulo cosmic censorship), a black hole would 
form in the spacetime being constructed. We note that, 
as is the case with the usual Schwarzschild coordinates 
for a spherically symmetric black hole, polar-areal coor- 
dinates cannot penetrate apparent horizons, and in fact 
become singular as they come "close to" black-hole re- 
gions of spacetime, where 2M{t,r)/r 1. This fact 
does not present a problem in the study of critical be- 
havior in our models, since the critical solutions per se 
have maxr {2M(t, r)/r} bounded away from 1. 



B. Regularity and Boundary Conditions 



In addition to the above equations of motion, appropri- 
ate regularity and boundary conditions are needed. At 
the origin, r — 0, regularity is enforced via 



(21) 
(22) 

(23) 
(24) 
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= O(r'), 


H(t,0) 


= O(r'), 


$(t,0) 


r 0(r' 
- \ 0(r) 



^) for I > 1, 
for I = 0. 



(25) 
(26) 



In the continuum, our equations of motion are to be 
solved as a pure Cauchy problem, on the domain t > 0, 
r > 0, with boundary conditions at spatial infinity given 
by asymptotic fiatness (i.e. that the matter fields vanish, 
and that the metric becomes that of Minkowski space- 
time, as r — > oo). Computationally, we solve an ap- 
proximation to this problem on a finite spatial domain 
< r < Tniax, where rmax is some arbitrary outer ra- 
dius chosen sufficiently large that we are confident that 
the numerical results do not depend significantly on its 
precise value. At the outer boundary, then, the following 
condition for a is imposed: 



a(i, Tinax) a(t,r^i,^) = 1. 



(27) 



This can be viewed as simply providing a convenient nor- 
malization for a, since given a solution, a, of the slicing 
equation (|18p . ka is also a solution, where k is an arbi- 
trary positive constant. We note that although we have 
used (|?7|) in order to perform the calculations, a different 
normalization convention — i.e. a different, and time de- 
pendent, choice of k — has been used in order to perform 
the analysis of the solutions. Specifically, in the analysis 
we have used central proper time T defined by: 



T = 



a{i,0)di. 



(28) 
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This definition of time has a natural geometrical interpre- 
tation since r = is invariantly defined by the symmetry 
of the spacetime. For the scalar field variables, 11 and 
$, approximate outgoing-radiation boundary conditions 
(Sommerfeld conditions) are used: 

f + f (t,w) + ii^ll^ . 0, (29) 

dt dr r,nax 

^(i, fmax) + ^(i, fmax) H = 0. (30) 



Family Form of initial data, V-'CO^ ^) 



(a) 



A exp(-(r-ro)Vcr^) 



A 



(b) 



-2A(r - ro)/g^exp (-(r - ro)Vg^) A 



(c) A (atan(r — Tq) — atan(r — tq — u)) A 



TABLE L Families of initial data and the parameter p that 
is tuned to generate a critical solution. 



An important point in the derivation of the equations 
of motion is the fact that the eigenfunctions in ^ are dis- 
crete and the allowable values of I are only non-negative 
integers. Once the equations are obtained we have re- 
laxed that constraint and have allowed / to take non- 
negative real values. The solutions corresponding to non- 
integer values of I would have some degree of irregular- 
ity at the origin depending on the particular value of I 
chosen. This implies that only some finite number of 
derivatives with respect to r will be defined at r = 0. In 
our particular numerical implementation, which assumes 
that second derivatives of the variables are defined, we 
have been able to study the evolution of these systems as 
long as 2 > 3. 



Initial Data (F) Family Parameters 



1 


(a) 


ro = 70.0, 


a = 


5.00 


2 


(b) 


ro = 70.0, 


a — 


5.00 


3 


(c) 


ro = 70.0, 


a = 


5.00 


4 


(a) 


ro = 40.0, 


a = 


10.0 


5 


(a) 


ro = 40.0, 


a — 


5.00 


6 


(a) 


ro = 70.0, 


a — 


10.0 



TABLE IL Initial data used in our investigations. The family 
labels are defined in Table I. 



B. Families of Initial Data 



III. RESULTS 

A. Numerics 

We solve equations pS)) , ([TB]) for the scalar field gradi- 
ents, equations p?)l . (fT8)) for the geometry, and use ([13]) 
to reconstruct the field t/i. The system is approximated 
using second order centered finite difference techniques, 
and coded using RNPL J5j. Numerical dissipation of 
the Kreiss-Oliger [l6l | variety was included to damp high 
frequency modes, and it should be noted that this par- 
ticular type of dissipation is added at sub-truncation er- 
ror order, so does not affect the overall accuracy of the 
scheme as the mesh spacing tends to 0. For the cur- 
rent computations, the damping terms were most useful 
in regularizing the truncation error estimation procedure 
that occurs when adaptive mesh refinement (AMR) tech- 
niques are used. It was also crucial to impose the correct 
leading-order regularity conditions close to the origin, 
r = (equations ([251) -(ESI)) 7 in order to keep the solution 
regular during the evolutions. Most of the calculations 
were done on a fixed uniform spatial grid r^ — [j — 1) Ar, 
j = 1, 2, • • • , J, J = 1 -I- rmax/Ar with a typical number 
of grid points J = 1025, and the outer boundary of the 
computational domain typically at rmax — 100. For small 
values of the angular momentum parameter — specifically 
for I < 2 — an AMR algorithm based on that described 
in [l3| was used. 



Our study involved the evolution of 6 different one pa- 
rameter families of initial data, each defined by an initial 
profile -0(0, r) as listed in Table U with specific values 
of the parameters appearing in the profile definitions as 
given in Table [TTl In addition to -0(0, r), we need to pro- 
vide n(0,r) to complete the specification of the initial 
data. In all cases we chose 11(0, r) to produce an approx- 
imately in-going pulse at the initial time: 

n(0,r)=<i>(0,r) = |^(0,r). (31) 

As previously mentioned, all of the initial data families 
listed in Table[T]have a single free parameter, p, and, as is 
the usual case in studies of black hole critical phenomena, 
for any given family we observe two different final states 
in the evolution, depending on the value of p. For values 
oi p > p* the maximum value of 2M{t, r)/r approaches 1 
implying that an apparent horizon is about to form. On 
the other hand if p < p* the scalar field completely dis- 
perses, and leaves (essentially) flat spacetime in its wake. 
The solution that arises as p — > then represents the 
threshold of black hole formation and, by definition, is 
the critical solution. We note that these critical solutions 
are not t — > 00 end-states of evolution; rather they persist 
for only a finite amount of time, and, in fact, are unstable, 
heuristically representing an infinitely fine-tuned balance 
between dispersal and gravitational collapse. 
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C. Analysis 

We have calculated p* for the different families of ini- 
tial data described above, and for different values of 
via bisection (binary search), tuning p in each case to a 
typical precision of {p — p*) jp k, 10~^^ (which is close to 
machine precision using 8-byte real floating point arith- 
metic). 

As in the case for I = (where the equations of mo- 
tion reduce to those for a single, non-interacting mass- 
less scalar field, as studied in ^17,]), the critical solutions 
for values of I < 9.5 are apparently discretely self sim- 
ilar (DSS). DSS spacetimes are scale-periodic, meaning 
that any non-dimensional quantity, Z, obeys the follow- 
ing equation for some specific values of the parameters 
A and T*: 



Z ((T -T*),r)=Z {e"^{T - T*), e"'^r) 



(32) 



where T is central proper time as defined by (|28p . and 
T* is the "accumulation time" of the self-similar solu- 
tion. In ((32|) the integer n denotes the "echo" number. 
We also note that due to the discrete ip —ip invari- 
ance that is exhibited both by the equations of motion 
as well as the critical solutions themselves, if A is the 
echoing exponent for which formula (j32p is satisfied with 
Z{T,r) = ip{T,r), then the geometric quantities a{T,r), 
a(T, r), 2M{T,r)/r obey ([32| with an echoing exponent 
A/2. 

In order to extract A from our calculations, we use the 
observation that certain geometric quantities will achieve 
(locally) extremal values on the spatial domain at dis- 
crete central proper times T„ given by 



Tn - T* = (To - T*) e 



TiA/2 



(33) 



where Tq is the time at which one starts counting 
the echoes. Specifically, A and T* have been com- 
puted by a least squares fit for the times r„ at which 
maxr {2M(t,r)/r} achieves a local maximum in time, 
i.e. by minimizing: 



N 
n=l 



D. Results 

Table IIIII summarizes the values of A; we have esti- 
mated using this procedure; the data are also graphed in 
Fig. [T] Again, note that the reported values for A; have 
been calculated using central proper time T instead of 
proper time at infinity (the parameterization used in the 
numerical evolutions per se). Also the reported uncer- 
tainties have been estimated from the deviations in the 
A; values computed across the the six different families 
of initial data. The first entry in Table IIIII (/ = 0) cor- 



/ 


A; 


7/ 





3.43 ± 0.05 


0.376 ± 0.003 


1 


0.460 ± 0.002 


0.119 ± 0.001 


2 


0.119 ± 0.003 


0.0453 ± 0.0002 


3 


0.039 ± 0.001 


0.020 ± 0.001 


3.5 


0.0224 ± 0.0009 


0.0127 ± 0.0008 


4 


0.0132 ± 0.0008 


0.0082 ± 0.0008 


4.5 


0.0077 ± 0.0007 


0.0052 ± 0.0006 


5 


0.0044 ± 0.0007 


0.0033 ± 0.0005 


5.5 


0.0026 ± 0.0006 


0.0020 ± 0.0005 


6 


0.0015 ± 0.0005 


0.0013 ± 0.0005 


6.5 


0.0009 ± 0.0005 


0.0008 ± 0.0005 


7 


0.0006 ± 0.0004 


- 


7.5 


0.0004 ± 0.0004 


- 


8 


0.0003 ± 0.0004 




8.5 


0.0002 ± 0.0003 




9 


0.0002 ± 0.0004 




9.5 


0.0002 ± 0.0003 





TABLE III: Summary of the properties of the critical solu- 
tions computed for different values of I. Note that both the 
echoing exponents. A;, and the mass scaling exponents, 7;, 
rapidly decrease as I increases. Quoted errors have been esti- 
mated from the variation in values computed across the differ- 
ent families of initial data. Values of Aj have been calculated 
using central proper time normalization of the lapse function, 
which is the natural normalization for type-II critical behav- 
ior. For I > 6.5 we have not been able to calculate 7; due to 
lack of numerical precision. Note that the I = data agree 
with the original values calculated in [TtII. and that the I = 1 
data agree with values calculated in [iSj and [l9| using models 
of completely different origin. 



responds to the original case studied in [131 ■ The second 
one (Z — 1) is apparently the same solution found for the 
self-gravitating collapse of an 5*0(3) non-linear a model, 
assuming a hedgehog ansatz [l^, Interestingly, the 
values for Ai and 71 also agree quite well with the val- 
ues obtained from the study of the axisymmetric collapse 
of a complex- valued scalar field with azimuthal quantum 
number m — 1 11], where values A « 0.42 and 7 » 0.11 
are quoted. However, in the model considered in ^1¥\, the 
overall solution is clearly different because it is not spher- 
ically symmetric. The remainder of the solutions (for the 
other values of /) are, to the best of our knowledge, new. 

Systems exhibiting type II critical behavior, where 
the critical solution is self-similar, generally also exhibit 
power-law scaling of dimcnsionful quantities in near- 
critical evolutions. For example, we can expect the black 
hole mass, Mbh, to scale as 



Af, 



BH 



Cip-p*) 



(35) 
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-1 
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I.D. 


1) 


• 


I.D. 


2) 


A 


I.D. 


3) 


• 


I.D. 


(4) 


■ 


I.D. 


(5) 




I.D. 


(6) 




FIG. 1: Values of logjQ (A;) versus I. In this figure we can see 
that A; decreases almost exponentially with I. The different 
lines represent different families of initial data. Assuming 
universality, the differences between the values calculated for 
the different families provides one measure of error in our 
determination of A;. 



FIG. 2: Values of logj^Q (7;) versus where 7; is the scaling 
exponent defined by (|35}. As for the case of the echoing 
exponent, Aj, 7; also decreases approximately exponentially 
with I. We note that due to lack of numerical precision we 
can only reliably compute 7; for / < 6.5 



for super-critical evolutions as p ^ p* Here C is a 
constant that depends on the family of initial data while 
7/ is a universal exponent for each value of I, i.e. indepen- 
dent of the specific initial data family used to generate 
the critical solution. We have observed such scaling in at 
least some of our computations, but, following Garfinkle 
and Duncan [20| have found it more convenient to extract 
7/ by monitoring the maximum value of the trace of the 
stress tensor, T, which, from the Einstein equations, is 
proportional to the maximum value of the Ricci curva- 
ture. On dimensional grounds T (defined by (H^) and R 
should both scale with an exponent —27. This technique 
has the advantage of being more precise than a strategy 
based directly on ((35)) since we can calculate the trace of 
the stress-energy more accurately than the mass of the 
black hole formed, and can perform the computation us- 
ing sub-critical evolutions, where the gradients of field 
variables generally do not become as large as those in 
the super-critical cases. The values of 7; as a function of 
I are listed in Table Hill and are plotted in Fig. [21 

As is characteristic of type-II critical solutions ex- 
hibiting discrete self-similarity, 2M(t,r)/r oscillates at 
higher frequencies and on smaller spatial scales during 



^ In accord with the results in Q and Q], we expect small ampli- 
tude oscillations with period A to be superimposed on the scaling 
law H35|l . We have, however, made no attempts to measure this 
effect in the current work. 



the course of an evolution in the critical regime. As has 
already been noted, as I increases, the echoing exponent 
A/ decreases rapidly. This can be observed in Fig. HI 
where the evolution of the maximum in r is shown as a 
function of time for four different values of L. 

In addition, also in Fig. |3l we observe that the max- 
imum and minimum values between which the spatial 
maximum of 2M{t,r)/r oscillates increase with I (this 
fact is shown for all values of I in Fig. U) indicating that 
the critical solutions are becoming increasingly relativis- 
tic as the angular momentum barrier becomes more pro- 
nounced. The amplitude of the oscillations between these 
extremal values decreases since min^ {2Af (i, r)/r} in- 
creases more rapidly than max^ {2M{t, r)/r} (see Fig.|4l). 

The assumption that the critical solutions arc indepen- 
dent of the initial family of initial data implies that the 
spatial profiles at the same moment during the oscilla- 
tion for two different families of initial data are the same 
up to some rescaling of the radial coordinate. In Fig. [51 
we show a check of the universality of the spatial profile 
for the solutions computed with 1 = 9. Specifically we 
compare the spatial profiles at times T„, times at which 
the local maximum in time is achieved during criticality, 
for different families F of initial data F = 1, 6 given in 
Table HIl In order to compare profiles we rescaled the ra- 
dial coordinate by a constant Kp, which depends on the 
family of initial data. These constants are chosen in such 
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FIG. 3: Evolution in time of the maximum in r of the function 
2M{t,r)/r for four different critical solutions with increasing 
value of Z (/ = 0, Z = 1, Z = 2 and 1 = 4). The plot shows the 
evolution during the period of time when each solution shows 
discrete self similarity. The time coordinate is rescaled by A; 
for visualization purposes and is shifted so that the function 
values coincide at t = 0. We note how the solutions tend to 
periodicity with increasing values of I. 
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FIG. 4; maxt {maxr {2Af (f, r)/r}} in the critical regime 
as a function of I (solid line) and the same for 
mint {maxr {2M(t, r)/r}} (dashed line). We see how both 
the maximum and minimum values of 2M/r increase with I. 
On the other hand the amplitude of oscillation, given by their 
difference, apparently tends to zero with increasing I. 



FIG. 5: In the top pane we show the spatial profiles (in the 
region of self-similarity) of the scalar field ■;/) for different fami- 
lies of initial data, but for fixed angular momentum parameter 
1 = 9. In particular we show the solutions calculated from 
initial data types F = 1, 6 (see Table |TI| at times when tpp 
reaches maximum amplitude. Each solution is shifted by an 
amount proportional to its family number for better visual- 
ization, with F = 1 the bottom curve, and F = 6 the top. 
The r coordinate is rescaled for each family by a constant 
factor Kf, which is family dependent, in such a way that the 
difference with respect to the profile obtained for the initial 
data labeled with F = 1 (for which we consider Ki = 1) is 
minimized. In the bottom pane we show the differences be- 
tween the rescaled profiles for F = 2, and the profile for 
F = 1, divided by the £2 norm of the solution. The maximum 
relative difi'erence is of the order of a few percent, providing 
strong evidence that the critical solution is universal. 



a way that the of the difference of the profiles 

with respect to the one with F — 1, which is considered 
to have Ki = 1, are minimized. We observed that the 
maximum of the relative difference, i.e. the difference di- 
vided by the £2-iiorm of the solution, is of the order of 
a few percent, providing strong evidence for universality. 
Similar differences have been observed for other values of 
the angular momentum parameter. 

Empirically, we have also found that, as we increase 
I within a family of initial data, although A; ^ and 
T;* 0, the product T;*A; appears to asymptote to a 
finite value. Note that ostensibly this product is family 
dependent (see Fig. (5]), but again that all DSS type-II 
critical solutions are universal only up to a global scale 



The ^2-norm of a vector u defined as = y ^'f/^- 



8 



30 



25 - 



20 



15 - 



10 



I.D. (1) 
•I.D. (2) 

U M.D. (3) 
•I.D. (-1) 
■I.D. (5) 

OLD. (6) 

* . f : : 



■ ^ - ^ - A - ^- - A - 4 - -A^ - A - 



FIG. 6: T;* A; as a function of L The fact that these products 



remain finite as T,* 



oo and A; 



is evidence that the 



critical solutions tend to a periodic solution in the limit 



transformation (r, t) — > (fcr, kt), with k an arbitrary pos- 
itive constant. Choosing k = k{l) for each of the famihes 
so that maxr {2AI{t, r)/r} is attained at some fiducial ra- 
dius ro, and considering the case / = 10, we find that the 
normahzed asymptotic oscillation frequency, /o, defined 
by 



that the maximal value is attained at times T„: 

r„ = - (T*A) n. (38) 

As shown in Figs. [7]and[8l from our calculations for I — 
10, we cannot ascertain whether the solution is discretely 
self-similar with A; very small (< 0.0002), or periodic 
with period r = r*A. 
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FIG. 7: Fit of the times T„ at which maxr {2M{t,r)/r} 
reaches its maximum in time (triangles, left scale) assuming 
a periodic ansatz. Initial data type F = 1 was used with 
angular momentum parameter I — 10. We also plot the resid- 
uals of each data point with respect to the best fit (pentagons, 
right scale). 

Naively at least, we expect that for ^ > 10, distin- 
guishing between discrete self-similarity and periodicity 
would become even more difficult. However, it is worth 



/o = ro/(r*A) =4.35 ±0.01 



(36) 



agrees for all families to better than 1%. Again, the 
quoted uncertainty is estimated from the variation of 
/o across the different families of initial data. We note 
that for I — 10 the near-critical solution stays at a near- 
constant radial position; our spatial resolution is insuffi- 
cient to resolve the small changes associated with the 
extremely small value of A;. The radial location of 
maxr {2M{t,r)/r} in this regime is the value of vq that 
we have used in 

We also note that the observation that /o is apparently 
well defined and unique (up to the usual rescalings asso- 
ciated with type-II critical solutions), is consistent with 
the empirical observation that as I increases, the critical 
solution becomes ever closer to a periodic solution. In 
particular, for a periodic solution we have A — » 0, and 
then 

Tn-T* = (To - T*) e"'^ « (Tq - T*) (1 + nA) 

« -{T*A)ji-T\ (37) 

where Tq represents the loosely defined time demarking 
the onset of the critical regime (and whose precise value 
is clearly irrelevant in the limit T* —>■ oo) which implies 



G 




5 10 15 
n 

FIG. 8: Fit of the times r„ at which maxr {2 M{t,r)/r} 
reaches its maximum in time (triangles, left scale) assuming 
a self-similar ansatz. As in the previous plot, initial data 
type F = 1 was used with angular momentum parameter 
/ = 10. Again, we also plot the residuals of each data point 
with respect to the best fit (pentagons, right scale). Notice 
that the errors in the fit are of the same order as the errors in 
the fit that assumes periodicity (Fig. [7]), indicating that from 
our numerical results we are unable to distinguish between 
the two types of solutions for I > 10. 
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noting that for Z = 20 we have not yet seen evidence for 
(almost)-periodicity, with period T*A, but have instead 
seen a more comphcated structure near criticaUty that is 
not yet understood. 



IV. CONCLUSIONS 

In this paper, we have discussed the results for a model 
that incorporates some of the effects of angular momen- 
tum in the context of critical gravitational collapse. A 
new family of spherically-symmetric critical solutions, 
(black hole threshold solutions) labelled by an angular 
momentum parameter, I, has been found. These solu- 
tions have similar properties to those for the I = case 
originally studied in specifically, the solutions ex- 

hibit discrete self-similarity, and have scaling laws for 
the values of dimensionful quantities in evolutions close 
to criticality. We have calculated the /-dependence of the 
echoing exponents A; , and the mass-scaling exponents ji , 
finding that both decrease rapidly with increasing Z, (at 
least up to Z « 10). Moreover, we have argued that as I 
increases, the critical solution approaches a periodic evo- 
lution. 

Together with the results of [iiIj our findings suggest 
that certain models of collapse may generically admit 
countable infinities of critical solutions, each member of 
which can be characterized by distinct near-origin regu- 
larity conditions (such as (|24j|26p l that are preserved by 
dynamical evolution. 

As we explained in the introduction, we expect that 
7; = l/Xi where A; is the Lyapunov exponent associated 
with the single unstable mode of the critical solution for 
angular momentum parameter I. Therefore since 7; ^ 
with increasing Z, we apparently have A; — > oo. This 
has the interpretation of increased stability of the critical 
solution for increasing Z, i.e. the period of time that a 
solution can remain close to criticality (for a fixed amount 
of fine tuning) increases with Z. We believe that this 
can be interpreted as an effect of the angular momentum 
barrier which (partially) stabilizes the collapse to black 
hole formation. 
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APPENDIX A 

We wish to show that the stress energy tensor T^'^'^f, — 
XI^TC")";, is independent of 6 and 0, where rC")''^ 
is the stress energy tensor associated with the solution 
ij^''\t,r)Qim{0,4>), with the same function -0^'^ (i, r), for 
each value of m. For the scalar field, the tensor Tab can 
be written in terms of the solutions to the wave equation, 
^, as 

T\ = 5"'*,c*b - ^<5%/"*,c*d , (Al) 

and if 

Y,M\ = g-(V-W(i,r)gi™(0,0))_^ 

m 

x(V^'^(i,r)g,™(0,0)),, (A2) 

is independent of 6*, 0, then so is T^'^^f,. 

Using the definition of the Qim this can be written in 
terms of the as 

x(^(')(t,r)i^„,(0,0)),. (A3) 
We can write this in terms of the Green's function 

p{9, 0, 9\ 0') = J2 mim{e', 4>') (A4) 

m 

in the limit as 6' ^ 9 and (/>' — > (j). 

In bra-ket notation, this is just the operator 

P = ^ |Zto) (Zto| (A5) 

which commutes with all of the angular momentum op- 
erators. 

[L„P] = 5][i.,|ZTO) (Zm|] 

— y^(m I Im) {lm\ — \ Im) {lm\ m) 

= 0, (A6) 

[L^ + iLy,P] = ^(V^(Z + l)-m2-m|Zm + l) (Zm| 

m 

-\lm) {{L,^iLy)\lm))^^ 
= (yijl + 1) - - TO |Zm -I- 1) (Zm| 

m 

+ 1) - m'^ + m\ Im) (Zto-1|) 
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(A7) and similarly the term 



= 0. 

Thus Y.m ^m(^' (l})Yim{e', (j)') must be a function of the 
only rotation invariant function of 0, (p, 0\<p\ which is the 
angle 8 defined by 

cos(e) = cos(6') cos(6l') + sin(6') sin(6l') cos(0 - (/>')-(A8) 

is the angle between the two unit vectors with direc- 
tions 0,0 and d' i^l respectively. Since P depends only 
on 8 we can choose 61 = to evaluate it, which gives 

= y;S(o,o)rio(e,o) 

= ^^'Kcos(8)). (A9) 

The various components of the tensor M are of three 
types: ones with no derivatives with respect to d or 
(\) (eg Mtt), those with one derivative, (for example 
Mtg) and those with two (eg, Mge)- The ones with 
no derivatives will be functions of lime'_0/_»e^0 P — 
\J {21 + l)/(47r) Pi{l) which is clearly independent of 9, (j). 
The terms with one 0, 4> derivative will be functions of 

hm^a9,0P(cos(8)) = hm^P/{sin(8),sin2(6')sin((?!)-0')} 
= 0, (AlO) 



Mg^ oc \im dgd^' Pi{cos{<d)) (X lim sin(0 — 0') 

— *0 (f) — *(j>' 

= 0. (All) 
Thus the only two terms remaining are 

Mge oc Van dede'Piicoiiie -d')) 

9-^9' 

= P/(cos(0))(-cos(0))(-l) 

= Pli^), (A12) 
M^^, oc P/(l)sin(0)2. (A13) 



Thus the non-zero components of the tensor M are 
Mtt, Mtr, Mrr, M^^ = SIT? {9)1^09, with Only M^^ hav- 
ing 9 dependence. Thus, M'^f, will be independent of 9, (j) 
and therefore so will T^'^^f,, as required. 

In addition, from the equations of motion for the indi- 
vidual fields '^im, each of the energy momentum tensors 
for given /, m is conserved in the overall spherically sym- 
metric spacetime. Thus, so is their sum over m for any 
given and we have 

rW-^M^O. (A14) 
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